!
!  photometry - photometry
!
!  Copyright © 2010 - 2011 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.


program photometry

  use ifits
  use mdaosky
  use mdaofind
  use mdaofotometr

  implicit none

  integer, parameter :: clen = 4096
  character(len=clen) :: file,line,record
  character(len=clen) :: output = ''
  real, dimension(12) :: apertures
  real, dimension(2) :: ring
  real :: lothresh = 7.0
  real :: threshold = 30.0
  real :: fwhm = 2.5
  real :: saturation = 6.4e4
  integer :: maxsky = 10000
  real :: shrplo = 0.2, shrphi = 1.0
  real :: rndlo = -1.0, rndhi = 1.0
  real :: readns = 1.0
  real :: phpadu = 1.0
  integer :: i,cc,eq,istat


  forall( i = 1:size(apertures) )
     apertures(i) = 2.0*exp(0.1567271184*((i-1)*1.570796327))
     ! logarithmic spiral, pitch 8.9073 deg
  end forall
  ring = (/20.0, 30.0/)


  do
!     read(*,*,end=20) sb

     read(*,'(a)',end=20,iostat=istat) record
     if( istat /= 0 ) then
        write(*,*) 'Unexpected string on input.'
        stop 21
     end if

     cc = index(record,'#')
     if( cc == 0 ) then
        line = record
     else
        line = record(1:cc-1)
     end if

     eq = index(line,'=')

     if( index(line,'FILE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) file
        if( istat /=0 ) then
           write(*,*) 'Invalid filename format.'
           stop 21
        end if

        write(*,*) trim(file),':'
           
        if( output == '' ) then
           write(*,*) 'An output file to store photometry not specified.'
           stop 21
        end if

        call finit(output,output==file,apertures,ring)
        call ffind(file,output,maxsky,saturation,fwhm,lothresh, &
             threshold,shrplo,shrphi,rndlo,rndhi,readns, phpadu)
        call faphot(file,output)

     else if( index(line,'OUTPUT') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) output
        if( istat /=0 ) then
           write(*,*) 'Invalid output filename format.'
           stop 21
        end if        

     else if( index(line,'FWHM') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) fwhm
        if( istat /=0 ) then
           write(*,*) 'Invalid FWHM format.'
           stop 21
        end if

     else if( index(line,'THRESHOLD') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) threshold
        if( istat /=0 ) then
           write(*,*) 'Invalid threshold format: ',trim(line)
           stop 21
        end if

     else if( index(line,'LOWER THRESHOLD') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) lothresh
        if( istat /=0 ) then
           write(*,*) 'Invalid lower threshold format.'
           stop 21
        end if

     else if( index(line,'ROUND LOWER') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) rndlo
        if( istat /=0 ) then
           write(*,*) 'Invalid lower round format.'
           stop 21
        end if

     else if( index(line,'ROUND HIGHER') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) rndhi
        if( istat /=0 ) then
           write(*,*) 'Invalid higher round format.'
           stop 21
        end if

     else if( index(line,'SHARP LOWER') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) shrplo
        if( istat /=0 ) then
           write(*,*) 'Invalid lower sharp format.'
           stop 21
        end if

     else if( index(line,'SHARP HIGHER') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) shrplo
        if( istat /=0 ) then
           write(*,*) 'Invalid higher sharp format.'
           stop 21
        end if

     else if( index(line,'SATURATION') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) saturation
        if( istat /=0 ) then
           write(*,*) 'Invalid saturation format.'
           stop 21
        end if

     else if( index(line,'PHPADU') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) phpadu
        if( istat /=0 ) then
           write(*,*) 'Invalid photons per ADU format.'
           stop 21
        end if

     else if( index(line,'READ NOISE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) readns
        if( istat /=0 ) then
           write(*,*) 'Invalid read noise format.'
           stop 21
        end if

     end if

  end do

20 continue

  stop 0

end program photometry
